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Abstract 

Stability margin for multiloop flight control systems has become a critical 
issue, especially in highly maneuverable aircraft designs where there are in- 
herent strong cross-couplings between the various feedback control loops. To 
cope with this issue, we have developed computer algorithms based on non- 
differentiable optimization theory. These algorithms have been developed for 
computing the Multivariable Stability Margin (MSM). The MSM of a dynam- 
ical system is the ”size” of the smallest structured perturbation in component 
dynamics that will destabilize the sytem. These algorithms have been coded 
and appear to be reliable. As illustrated by examples, they provide the basis 
for evaluating the robustness and performance of flight control systems. 


1 Introduction 

Accurate knowledge of the dynamical model associated with the design of modern 
flight control system is becoming more difficult to obtain. This is especially true for 
the design of the next generation fighters where many of the performance specifica- 
tions go beyond the capability of the aircraft currently in service. Robust control 
analysis methods have received considerable attention in recent years as a possible 
solution to the problem of controlling systems for which the given model contains 
significant uncertainty [Saf 1] [Doyle2]. The central feature of these methods is their 
effectiveness in handling an unknown-but-bounded class of plants, instead of the nom- 
inal plant only. 
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The uncertainty in the nominal plant arises from several different sources: For 
gain-scheduled aerospace vehicle control systems, typical uncertainties in the plant at 
each design point consist mainly of modeling errors due to uncertain aerodynamic co- 
efficients, linearization, model reduction, neglected dynamics, time-delays, etc. Aero- 
dynamic coefficients developed from wind tunnel testing or computational fluid dy- 
namics usually are different from those obtained from actual flight data. Linearization 
will also affect the nominal plant behavior. Nonlinear effects such as actuator satura- 
tion and rate limits are neglected altogether when a model is linearized. Parameter 
drift will also affect the nominal plant. There may be dynamical modes which are 
intentionally or unknowingly neglected. In addition phase loss which results from 
time delay also leads to an uncertainty bound. 

The uncertainty may be loosely classified as falling into two categories, structured 
and unstructured. Structured uncertainty arises from specific component or parame- 
ter variations. Two examples of structured uncertainty are variations in weight and 
drifting aerodynamic parameters. Unstructured uncertainty is any other sort of uncer- 
tainty which can be regarded as a frequency-dependent norm bounded perturbation 
matrix. High frequency modeling errors are one type of unstructured uncertainty. 
The linearization of the nonlinear equations of motion contribute to both classes of 
uncertainty. Actuator rate and position limits have a distinctive signature which can 
easily be isolated, so that they fall into the class of structured uncertainties. On the 
other hand, the effects of nonlinear kinematic terms can only be bounded, therefore 
necessitating an unstructured uncertainty representation. 


2 Robustness Measure 

Safonov [Saf 2], who built upon different, but related, conic-sector nonlinear stability 
theory work of Zames [Zames], reinterpreted the conic-sector stability concepts in 
order to deal with uncertainty and robustness issues. Safonov, Doyle and Fan, to 
name a few, have contributed to the continuing development in this area [Saf 3] 
[Saf 4] [Doyle2] [Fan]. 

Basically the robustness measure is done by lumping uncertain deviations from a 
nominal system M(s) into an uncertain matrix A(a) resulting in an uncertain feedback 
system with loop transfer function A (s)M(s) as shown in figure 1. 

Then, the Multivariable Stability Margin (MSM), ”K m n } is defined as the smallest 
stable, norm-bounded perturbation A(s) that can destabilize the system. While K m 
is in general difficult to compute, a reasonably tight lower bound K m theoretically 
can be computed using diagonally scaled singular values [Saf 3] [Doyle2] [Fan]. The 
plot of K. m (w) vs. frequency identifies tolerable levels of parameter uncertainty as a 
function of frequency. 
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Figure 1: Robustness analysis model. 

For unstructured uncertainty, the maximum singular value has been shown to be 
useful in bounding the multivariable stability margin. However, the bound can be 
very conservative in the case of structured uncertainty. The singular value analysis 
will attempt to find the worst direction of the uncertainty that in reality impossible to 
exist. To deal with the case of diagonally structured A, Safonov [Saf 4] introduced the 
two-sided structured Multivariable Stability Margin (MSM), denoted K m , and Doyle 
[Doyle2] introduced the term Structured Singular Value (SSV), denoted fi, to describe 
the reciprocal, /t(A/(a)) = l/if m (Af(a)). Diagonal perturbations are quite general 
and flexible if one considers parametric uncertainties (e.g. aerodynamic coefficients). 
Traditionally one defines K m and fi for ’’two-sided” magnitude-bounded uncertainties 
which may be either positive or negative; but in cases where the sign of the uncertain 
Aj is known a priori one may modify the definitions of K m and p. accordingly. 

When the uncertainties are known to cover both positive and negative pertur- 
bations, the SSV of Doyle and MSM of Safonov provide a ’’tight (to within 15%) 
condition for robust stability. This condition is measured by representing directly the 
individual sources of uncertainties in the form of block diagonal perturbations. 

Definition: Given transfer functions G(a),a(s) and 6(s), we write 

G(a) € sector[a, 6] 


if 


|G(ju;) - C{jw)\ < |r( < ;’u;)| Vw 


where 


and 


C(s) = (a(a) + b(s))/2 

r(s) = (a(s)-b(s ))/ 2 


Assuming A/(a) and A(s) to be stable then the one-sided MSM, ” K mi n and, the 
two-sided MSM, n K ma n are defined by the following (see figure 2 and 3): 
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Figure 2: One-sided K m 
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Figure 3: Two-sided K m 


• One-sided K m : 

The system is stable for all A with 

Ai G sector[0, K mi ] V» = 1, (1) 


• Two-sided K m : 

The system is stable for all A with 

A i £ 3ector[-K ma ,K m ,} Vt = 1, •••,«. (2) 

For any diagonal matrix D, a practical upper bound on fi = 1/K mj is cr max (DM D~ l ). 
Further, it is known that for 3 or fewer A; ’s that the minimum over D of this 
upper bound is actually equal to fi [Doyle2]. Safonov and Doyle proved that the 
minimization problem of o’mo X (DM D' 1 ) is convex in D' = log(D), so that every local 
minimum is a global minimum. Furthermore, computational experience has shown 
that minimum of a max (DMD~ 1 ) over D is within 15% of fi. So, we choose to work 
with this upper bound, to calculate the reciprocal of two-sided MSM. Practical 
upper bound for the reciprocal of one-sided MSM can be easily derived by using conic 
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Figure 4: K mr 


sector property of Zames [Zamesj; viz. K^] is bounded above by 

iC x = max[^minr>\ mat (DM D' 1 + (Z)MD -1 )*),0] (3) 

By modifying the former equation to include the permutation matrix fa, a less 
conservative bound for the two-sided real MSM, K mr , is given by 

£Lm r = rnax^max^minD^matiDMfaD -1 + (DMfaD' 1 )*), 0] (4) 

Here 

fa€$, * = 1, • • • ,2 n 

and $ is the set of all permutation nx n diagonal matrices. 

$ = diap[(±l), ,(±1)] 

The bound (4) is similar to the one proposed by Jones [Jones]: 

\min D max^\ max {DMfaD- 1 + ( DMfaD~ l y ) (5) 

Although, equation (5) will lead to more conservative bounds. Geometrically equation 
(4) is shown in figure 4. 

It should be mentioned at this point that several software packages are available 
to compute the two-sided MSM , however , they are not accessible to the authors to 
be evaluated. 

3 K m Computation 

The computation of <r max is straight forward using available numerical software (Lin- 
pack). For both one-sided, real and two-sided cases, we have to 6olve an optimization 
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problem. However for all of these problems, the analytical gradient is available, so 
accurate solution can be obtained. However, when several eigenvalues or singular 
values coalesce (i.e., have multiplicity greater than one) the function iB nondifferen- 
tiable (’’creases” produce direction-dependent derivatives), so that a more complex 
algorithm of computing a descent direction is required. Before actually solving either 
case, one can approximately prescale the system matrix M by substituting for M the 
matrix DMD -1 where D minimizes the Frobenius norm of DMD~ l [Osborne]. 

The monotonic transformation of P — ♦ D', with D = Exp(D'), transforms the 
problem into a well behaved convex optimization [Saf 5]. The initial guess for D was 
taken to be equal to the identity. This initial guess was used only for the first frequency 
value in the given range. The solution obtained for a particular frequency point was 
then used as an initial guess for the next value. Suppose that the largest eigenvalue 
is simple, then a descent direction is calculated directly using the Davidon-Fletcher- 
Powell technique. In the case that the largest eigenvalue has multiplicity greater than 
1 and the function iB not continuously differentiable, a generalized gradient is used 
to determine a descent direction. Once this is done, the minimal point can be found 
in the specified descent direction by using a well known ’’binary search” algorithm 
of Bolzano. These steps are repeated until the global minimum is located (i.e. gra- 
dient is zero). Convexity of cr^ oe (e £> 'Me _I> ') and X mam \{e D ' Me~ D ' -I- (e D> Me~ D ') m ) 
ensures that this procedure is convergent to the global minimum. These steps can be 
summarized as follows: 

1. Initialize Mi = M\ D[ = 0; k = 1. 

2. Scale M k+i = e £> *M*e -I> *; set D k+1 = 0 

3. Find the search direction 

- Davidon Fletcher Powell (DFP) deflected gradient. 

— DFP generalized gradient for multiplicity > 2 

4. Unidirectional search. 

- Method of Bolzano (Fig. 5). 

5. D' k+1 «— D' k+l + stepsize * search direction. 

6. k = k - (- 1; go to step 2. 

Step 5. of the algorithm involves varying the diagonal scaling matrix D' along a line 
by adjusting the scalar parameter stepsize. The size of e 0 *+» could approach the 
value of oo. To prevent this, as the stepsize grows, M k +\ , is repeatedly updated to 
Mk+\*~ e £>k + 1 A/e -£>t +» and stepsize and D' k+1 are reset to zero. This is done as often 
as needed to prevent numerical overflow when e £>t +> is evaluated. 
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Distance. 

Figure 5: Method of Bolzano. 


3.1 Generalized Gradient. 

This discussion is not at all self contained and only key results will be stated. An 
excellent reference for this algorithm is [Polak]. In the case where the greatest eigen- 
value/singular value has multiplicity greater than one, the function ceases to be dif- 
ferentiable. In this case the gradient is not defined and more complicated ’’generalized 
gradient” methods must be used to compute the descent direction. The generalized 
gradient at a nondifferentiable point is defined as the nearest point to the origin 
in the convex-hull of the set of directional derivatives at neighboring points; thus 
the computation of the generalized gradient at any point is itself a convex nonlinear 
programming problem. We employ an algorithm similar to that of [Doyk2, Polakl] 
to compute the generalized gradients of o’^ uut (e D Me~ D ) and A ma *|(e Me + 
(e D 'Me- D ')‘). Geometrically the algorithm is shown in figure 6 and summarized as 
follows: 

• Generalized Gradient is defined as: 

£ jv,(Co{v(*)| ||.|| = 1}) 


where 

Co(.) - the convex hull of the set (.). 

Nr(.) - the nearest point to the origin of the set (.). 
{V(*)| ||*|| = 1} - the set of directional derivatives. 

• Iterative algorithm for computing V flen : 

1. Initialize k = 1. 

2. Guess Xk and set yk = V(**). 
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Figure 6: Generalized gradient. 

3. Find xu+i by minimizing (yjj V(*fc+i)) subject to ||*fc+i|| = 1. 

4. Find y h+1 = Nr( Co(y fc , V(* fc+1 )). 

5. Increment k «— k + 1, go to step 3. 

3.2 Davidon-Fletcher-Powell Scaling. 

The unmodified generalized gradient determines a steepest descent direction. The 
steepest descent direction is simply minus the generalized gradient. Steepest descent 
usually works quite well during early stages of the optimization process but if the 
Hessian (second derivative) matrix has a large condition number, the method usually 
behaves poorly, and small zig-zagging steps, called ’’stitching”, take place (see fig. 
7). Stitching problems also occur when the multiplicity of <r max or X max is 3 or 
more. Therefore we use the Davidon-Fletcher-Powell (DFP) method to modify the 
generalized gradient in order to handle this phenomenon. This technique uses the 
previously calculated generalized gradient to estimate the Hessian and effectively 
rescale the function to make its Hessian better conditioned. This quadratic fit method 
requires fewer gradient evaluations and tends to converge faster. It should also be 
noted that the likelihood of stitching-induced premature termination of the algorithm 
(as can occur in the unsealed steepest descent technique) can be greatly reduced with 
the DFP scaling. 

4 Lateral Directional Flight Control Example 

4.1 Example 1 
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Figure 7: Minimization paths. 
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Figure 8: Axis systems and sign convention 
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Figure 9: Two-sided actuator uncertainty model. 
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Figure 10: One-sided actuator uncertainty model. 

Nonlinear elements of actuators can be treated as linear conic-sector elements with 
structured uncertainty. This uncertainty can be modeled as one-sided or two-sided 
uncertain gains within the actuator model. This can be shown through an example. 
Consider a saturation curve below. If the input size is always less than C , then the 
saturation element is equivalent to a gain element with a magnitude of one, however, 
if u exceeds C, one may model the saturation element as a two-sided uncertain gain 
A in parallel with a nominal gain of one as shown in figure 9. A better approach is to 
model the saturation element by a gain with a nominal value of one and a one-sided 
negative uncertainty as shown in figure 10. Clearly, the one-sided model will produce 
less conservative margins than the two-sided model. 

A design example is presented below in which MSM algorithm is asked to check the 
robustness of a typical lateral/directional flight control systems with respect to the 
actuator uncertainty (e.g. position saturation) and the reduction in the effectiveness 
of all control surfaces. The state-space matrices are given in figure 11. The controller 
uses roll rate, P, yaw rate, R , and the lateral acceleration, N y , for feedback (see figure 
12). By putting “extender wires” on the uncertainty blocks Aj and pulling them out 
into a separate “block”, one can check the system robustness. 

The plot of K mi , K. m3 and a max are in figure 13. Note that the cr moe and K m] 
which have the minimum values of .015 and .42 respectively are equal or less than 
K_ mi for all frequency, and therefore shown to be more conservative than one-sided 
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-3.4d-1 -4.8d-l -2.8d+2 3.2d+1 0.0d+0 
5.2d-3 -3.0d+0 4.6d-l -8.3d-6 0.0d*0 
1.5d-2 — 1 ,4d— 1 -1.9d + 0 -7.5d-7 0.0d*0 
0.0d+0 1.0d»0 0.0d*0 0.0d+0 0.0d+0 
0.0d+0 0.0d*0 l.0d*0 0.0d+0 0.0d+0 

-9.2d-4 -6.8d»* 

t .2d* 1 -7.2d-2 

3.8d-l 3.2d*t 

0.0d+0 0.0d*8 

O.Od+O 0.0d+8 

0.0d+0 5.7d+1 0.0d+0 0.0d+0 0.0d*0 

O.Od+0 0.0d+0 5.7d*1 0.0d+0 0.0d+0 

-6.3d-3 -3.6d-2 7.5d-1 5.1 d-8 0.0d+0 

o.od+o o.od+t 

0.0d*0 0.6d*0 

6.2d-2 1.H-I 



Figure 11: State-space matrices. 


u , 

Rellstk Compensation. 


$♦1.5 


ftpedel 


FL 

Compensation. 


A1 





A2 





A3 





Ad 


Rctuetor 


1/57.3 


.04S+1 


N * 

.645+2.5 

') * 

S 



>•- : 


i i 

♦ U 

-I — -d> 


actuator 

i 

t 

w 

1/57.3 

U 


.04S+1 




.63 


. 1 5S ♦ 1 


FB 


.49S(S+1K2(2S+t) 


<25*1X5+1) 


1.558 

.!$♦! 


N» 


•ircroft 

dynamics! 


R-ftX+BO 

T-CR+DD 


r» 

Compensation 


Figure 12: Lateral directional flight control with uncertainties at the controller output 
and plant input. 
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Figure 13: MSM at the plant input and controller output. 


structured stability margin K. mi . This results from the one-sided structured mul- 
tiplicative uncertainty that is not accounted for in the computation (i.e. nonlinear 
elements of actuators and gain reduction tolerance at the controller outputs). To 
properly account for the sign of the uncertainty and its structural information, the 
one-sided MSM was computed and it is shown to have a better robustness measure. 
2C m , has the minimum value of 0.677, indicating that the system can simultaneously 
tolerate at least a 67.7 percent reduction in the effectiveness of all control surfaces and 
the actuator inputs up to at least three times the saturation value without instability. 


4.2 Example 2 

The MSM’s minimal value K. mreak also can be used to quantify a control systems 
tolerance of simultaneous gain and phase variations at all the plant inputs and out- 
puts. This is done so the system has good stability robustness with respect to the 
uncertainties at the two actuator commands and the three sensor outputs (see figure 
14). These uncertainties come from various sources. Model accuracy deteriorates at 
higher frequencies due to unmodeled aeroservoelastic effect. Several potential error 
sources exist within the assumed perfect sensors. Model reduction of the actuators 
can also be considered as one of the effects at the plant inputs. 

Shown in figure 15 is the Bode plot of the MSM vs. frequency. The minimal 
value, denoted K m2ft k = .4196, gives an indication of the minimal size of structured 
perturbations required to destabilize the system or equivalently diagonal perturbation 
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Figure 15: MSM at the plant input and output. 
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as large as 41.96 percent can be tolerated at any frequency; at higher frequencies, 
perturbation magnitude as large as w/10 can be tolerated. 

5 Conclusion 

Computer algorithms for determining the multivariable stability margin u K m n have 
been developed. The algorithms provide a reliable tool for evaluating the robust- 
ness of control systems with significant gain and/or parameter uncertainties. The 
computation for the one-sided and two-sided structured stability margin were done 
using nondifferentiable optimization theory. Robustness analysis was performed on a 
typical lateral directional flight control system problem with large uncertainty. 
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